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ABSTRACT 

We conduct a survey of numerical simulations to probe the structure and appear- 
ance of non-radiative black hole accretion flows like the supermassive black hole at the 
Galactic centre. We find a generic set of solutions, and make specific predictions for 
currently feasible rotation measure (RM) observations, which are accessible to current 
instruments including the EVLA, GMRT and ALMA. The slow time variability of the 
RM is a key quantitative signature of this accretion flow. The time variability of RM 
can be used to quantitatively measure the nature of the accretion flow, and to differ- 
entiate models. Sensitive measurements of RM can be achieved using RM synthesis or 
using pulsars. 

Our energy conserving ideal magneto-hydrodynamical simulations, which achieve 
high dynamical range by means of a deformed-mesh algorithm, stretch from several 
Bondi radii to about one thousandth of that radius, and continue for tens of Bondi 
times. Magnetized flows which lack outward convection possess density slopes around 
-1, almost independent of physical parameters, and are more consistent with observa- 
tional constraints than are strongly convective flows We observe no tendency for the 
flows to become rotationally supported in their centres, or to develop steady outflow. 

We support these conclusions with formulae which encapsulate our findings in 
terms of physical and numerical parameters. We discuss the relation of these solutions 
to other approaches. The main potential uncertainties are the validity of ideal MHD 
and the absence of a fully relativistic inner boundary condition. The RM variability 
predictions are testable with current and future telescopes. 
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1 INTRODUCTION 

The radio source Sgr A* at the Galactic centre (GC) is 
now accepted to be a supermassi ve black hole (Mbh — 
4.3 x 1O 6 M : iGillessen et alJl200Sh . accreting hot gas from 
its environment [n e — 130cm -3 , ksT ~ 2keV at 1 arc 
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second: iBaganoff et~alll2003h . Interest in the Sgr A* accre- 
tion flow is stimulated by its remarkably low luminosity; 
by its similarity to other low-luminosity AGN; by circum- 
stantial evidence for pas t episodes of bright X-ray emission 
(|Revnivtsev et all |2004 but see lYusef-Zadeh et all l2007t ) 
and nearby star formation l|Levin fc Belo borodov 2003J) ; and 
foremost, by its status as an outstanding physical puzzle. 

Supermassive black holes are enigmatic in many re- 
spects; for the GC black hole (GCBH) the enigma is sharp- 
ened by a wealth of observational constraints, which per- 
mit detailed, sensitive and spatially resolved studies of 
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its accretion dynamics. Within a naive model such as 
Bondi flow, matter would flow inward at the dynami- 
cal rate from its gravitational sphere of influence, which 
at ~ 1" is resolved by Chandra. Converted to radiation 
with an efficiency nc 2 , the resulting luminosity would ex- 
ceed what is actually observed by a factor ~ 10 5 (T7/O.I) . 
This wide discrepancy between expectation and observation 
has stimulated numerous theoretica l explanations, includ- 
ing co nvection dNaravan et alJ 120001: iQuataert fc Gruzinovl 
l2000bh . outflow (|Blandford fc Begelmanl 1 19991) . domina- 
tion by individual stars' wind s ( Loebl |2004|). and con- 
duction iTranaka fc Menodl200d: I Johnson fc Quataertll2007l : 
ISharma et al.ll2008l ; IShcherbakov fe Baganofj|201Ch . 



1.1 Constraining the accretion flow 

Because many of its parameters are uncertain, the cen- 
tral density and accretion rate of the GCBH flow 
are n ot strongly constrained b y the its emission spec- 
trum (|Quataert fc Naravanlll999T ); the most stringent con- 
straint s come from observations of the rotation mea- 
sure (iQuataert fc Gruzinovl IgOOOah. now known to be 
roughly -5.4 x 10 5 radm _1 (|Marrone et all 1200?!) . Inter- 
preting this as arising within a quasi-spherical flow with 
magnetic fields in rough equipartition with gas pressure, 
and adopting the typical assumption that magnetic fields 
do not reverse rapidly, we derive a gas density nn ~ 
10 5 ' 5 cm" 3 (Rs / Rrd) 1 ^ 2 at the radius 7? re i which dominates 
the RM integral, namely where electrons become relativis- 
tic; see [JX] for more detail. If this radius is about 10 2 
Schwarzschild radii ( 1 2 _Rg ), as in the spectral models of 
IQuataert fc Naravanl (| 19991 ). then a comparison between 
this density and conditions at the Bondi radius Rb — 
0.053 pc indicates a density power law p oc r~ k with k = 
1.1 — 1.3;the derived value is rather insensitive to the black 
hole mass, the degree of equipartition, and the precise radius 
at which electrons become relativistic. (If rapid conduction 
causes electrons to be nonrelativistic at all radii, the implied 
slope is falls to 0.8.) 

An independent but weak constraint on k comes from 
recent multi-wave length observations of fla res in the emis- 
sion from Sgr A*. lYusef-Zadeh et al.l (2009) favor an inter- 
pretation in which these flares originate within regions in 
which electrons have been transiently heated and acceler- 
ated; using equipartition arguments they estimate a mag- 
netic field strength B ~ 13 — 15 G at 4 — 10 Schwarzschild 
radii, implying a total pressure P > 20dyncm~ 2 at those 
radii. Because P oc r~( h+1 \ a comparison to the conditions 
at Rb requires k > 0.6 — 0.8. This constraint could be vio- 
lated if the emitting regions were sufficiently over-pressured 
relative to the surround ing gas; however the sub sonic rate 
of expansion inferred bv lYusef-Zadeh et all l|2009h suggests 
this is not the case. 

The density power law k is an important diagnos- 
tic, both because it allows one to estimate the mass ac- 
cretion rate onto the black hole, and because k takes 
definit e valu es within proposed classes of accretion flows. 
iBondil (|l952h a ccretion and ADAFs (advection-dominated 
accretion flows, iNaravan fc~YHll994l ). in which gas under- 
goes a modifie d free fall, imply k = 3/2 and have long 
been ruled out (lAgo]||2000l ) by limits on the rotation mea- 
sure l|Bower et al.lll999l ~ CDAFs (convection-dominated ac- 



cretion flows, INaravan et alJ 120001 : IQuataert fc Gruzinovl 
2000b) and r elated flows like CDBFs (convect ion-dominated 
Bondi flows, llgumenshchev fc Nar avan 2002), in which con- 
vection carries a finite outward luminosity, all have k — 1/2 
outside so me small radius : otherwise, convection becomes 
supersonic l|Gruzinovll200ll ). 

Three classes of flows are known to have intermediate 
values, 1/2 < k < 3/2, as suggested by the observatons. 
One of these is th e ADIOS (advection-domina ted inflow- 
outflow solutions, iBlandford fc Begelmanl 1 19991 ). in which 
mass is lost via a wind from all radii within a rotating 
ADAF; however these flows appear to require that low an- 
gular momentum material has been removed from the axis. 
Another is a class of conductive flows, in which heat is car- 
ried outward by electrons a nd stifles accretion at large radii 
l| Johnson fc Quataerj|2007l 'l. A third consists of flows which 
lack any significant outw ard convective or conductive lu- 
minosity (|Gruzinovll200ll ). but are nevertheless hydrostatic 
rather than infalling; this behavior is seen within some nu- 
merical si mulations in which magnetiz ed g as is accreted, suc h 
as those o f llgumenshchev et al.l (20031 ) and lPen et all (|2003h , 
who termed the flow "magnetically- frustrated convection" . 

We are concerned with the last flow class, as it is phys- 
ically simple, realizable within simulations, and consistent 
with observational constraints. Whether it is physically rel- 
evant depends on the strength of conduction in the accretion 
flow, a question we return to in §[5] Although it is of interest, 
previous simulations do not suffice to make any quantita- 
tive comparisons between i t and the Sgr A* accretion flow. 
llgumenshchev et al.l (|2003h have already discussed several 
shortcomings which afflicted prior numerical work, such as 
(1) a lack of energy conservation during magnetic recon- 
nection and (2) simulation durations too short to capture 
steady states or secular trends. There are a number of other 
roadblocks: (3) Dynamical range: Rb is 10 5 Schwarzschild 
radii, but the largest simulations yet done have only a fac- 
tor of ~ 10 2 separating their inner and outer boundaries; 
(4) Resolution: numerical solutions are rarely close enough 
to the continuum limit to allow turbulent phenomena to be 
predicted with confidence; (5) Outer boundary conditions: 
although matter is presumably fed into the accretion flow 
by st ellar winds from the nuclear star cluster (jGenzel et al.l 
2003), the flow structure and magnetization of this gas is not 
well constrained; (6) Inner boundary conditions: the hole in- 
teracts with the flow in a manner which is not fully char- 
acterized, and which is likely to dominate the energetics; 
(7) Mass injection: stars within Rb produce fresh wind ma- 
terial, which have the potential to affect the final solution 
(|Loebll20"04l ): and (8) Plasma physics: close to be black hole, 
the flow is only weakly collisional, leading to effects such as 
anisotropic pressure and conduction, which may alter the 
natur e of fluid instabilitie s and the character of heat trans- 
port (Sha rma et al.l I2OO8I ). Potential deviations from ideal 
MHD become stronger as one approaches the event horizon, 
and are discussed further in section [5] 

In this paper we describe a numerical parameter sur- 
vey designed to partially overcome difficulties (l)-(5) in the 
above list, while making an educated guess regarding (7) and 
leaving (6) and (8) to future work. Specifically, we conduct 
three dimensional, explicitly energy conserving simulations 
to the point of saturation - often tens of dynamical times at 
Rb ■ We vary the dynamical range and resolution in order to 
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gather information about the astrophysical limits of these 
parameters, although they lie beyond our numerical reach. 
We push numerical outer boundaries far enough from Rb 
to minimize their effect on the flow, and we vary the con- 
ditions exterior to Rb in order to gauge the importance of 
magnetization and rotation in the exterior fluid. Our simu- 
lations obey ideal MHD, but are viscous and resistive on the 
grid scale for numerical reasons; we make no attempt to cap- 
ture non-ideal plasma effects. We do not account for stellar 
mass injection within the simulation volume. Our gravity is 
purely Newtonian, and at its base we have a region of accre- 
tion and reconnection rather than a black hole (although we 
are currently pursuing relativistic simulations to overcome 
this limitation). Our numerical approach is described more 
thoroughly below. 

By varying the conditions of gas outside Rb and by 
varying the allocation of grid zones within Rb we are able to 
disentangle, to some degree, physical and numerical factors 
within our results. We also compute integrated quantities 
related to the value and time evolution of RM, and draw 
conclusions regarding the importance of RM(f) as a powerful 
discriminant between physical models. 

We reiterate that our simulations have two simplifica- 
tions which could substantially change the behaviour. 1. Our 
black hole boundary condition is Newtonian. Since the deep- 
est potential dominates the dynamics and energy of the flow, 
a change in this assumption might alter the solution. 2. We 
assume ideal MHD to hold. As one approaches the black 
hole, the Coloumb collision rate is insufficient to guaran- 
tee LTE. Plasmas can thermalize through other plasma pro- 
cesses, but if these fail, strong non-ideal effects could dom- 
inate and lead to rapid conduction. These effects are both 
strongest at small radii, potentially modifying the extrapo- 
lation to the actual physical parameters. We address these 
issues in more detail in section \5\ 



From the above dimensional quantities we define several di- 
mensionless physical parameters. 

The adiabatic index is 7 = 5/3; the initial plasma-/? 
parameter, or ratio of gas to magnetic pressure, is 

8ttjpqc 2 s0 

po = — H2 — ; ( 2 ) 

we consider models with /3q = (1, fO, 100, 1000, 00) to cap- 
ture a wide range of plausible magnetizations. In our main 
sequence of simulations we adopt a uniform magnetic field 
BoS 

The initial velocity field is vo = (jo x f ) /r , where r is the 
separation from the black hole. The specific vector angular 
momentum is thus jo at the rotational equator, with solid- 
body rotation on spherical shells away from the equator. A 
dimensionless rotation parameter is therefore 

Rk _ ( joCs\ 2 ,„\ 
Rb ~ \Gm) ' (6) 

here Rk = jo / (GM) is the Keplerian circularization radius 
of the equatorial inflow. (Our flows never do circularize at 
Rk, both because angular momentum transport alters the 
distribution of j, and because gas pressure can never be 
neglected.) 

We impose mass accretion and magnetic field recon- 
nection within a zone of characteristic radius Ri n , de- 
scribed below, which introduces the dynamic range param- 
eter Rs/Rin- Because it sets the separation between small 
and large scales and the maximum depth of the potential 
well, this ratio has a strong influence on flow properties. 
One of our goals is to test how well the flow quantities at 
high dynamic range can be predicted from simulations done 
at lower dynamic range, as the dynamic range appropriate 
to Sgr A* is beyond what we can simulate. 



2 SIMULATION 

2.1 Physical setup and dimensionless physical 
parameters 

We wish our simulations to be reasonably realistic with re- 
gard to the material which accretes onto the black hole, but 
also easily described by a few physical and numerical param- 
eters. We therefore do not treat the propagation and shock- 
ing of individual stellar winds or turbulent motions, but take 
the external medium to be initially of constant density po 
and adiabatic sound speed c s o, and imbued with a character- 
istic magnetic field Bo and characteristic rotational angular 
momentum jo (but no other initial velocity). A Keplerian 
gravity field — GM/r 2 accelerates material toward a central 
"black hole" of mass M surrounded by a central accretion 
zone. The Bondi accretion radius is therefore 



We adopt the Bondi time is = Rb/c s o as °ur basic time 
unit; this is 100 years for the adopted conditions at Sgr A*. 
All of the initial flow quantities will evolve as a result of this 
during the course of the simulation, and we run for many 
Bondi times in order to allow the accretion flow to settle 
into a final state quite different from our initial conditions. 



2.2 Grid setup and numerical parameters 

We employ a fixed, variable-spacing Cartesian mesh in which 
the grid spacing increases with distance away from the black 
hole. To simplify our boundary conditions, we hold the spac- 
ing fixed within the inner accretion zone and near the outer 
boundary. The total box size is 4000 3 in units of the mini- 
mum grid spacing; however this is achieved within a numeri- 
cal grid of only 300 3 to 600 3 zones. Our grid geometry allows 
for a large number of long-duration runs to be performed at 
respectable values of the dynamic range, while avoiding co- 
ordinate singularities and resolution boundaries. These ad- 
vantages come at the cost of introducing an anisotropy into 
the grid resolution; however we have tested the code for con- 
servation of angular momentum and preservation of magne- 
tosonic waves, and found it to be comparable in accuracy to 
fixed-grid codes with the same resolution. Our grid expan- 
sion factor s = Sdxi/dxi takes one value for Xi < Rb and 
another, larger value for Xi > Rb', this allows us to devote 

1 We also investigated scenarios with Gaussian random field 
components, in which the dominant wavelengths were some mul- 
tiple of Rb; however we abandoned these, as such fields decay on 
a AlfVen crossing time, confounding our attempts to quantify the 
accretion flow, and we did not wish to add a turbulent driver to 
maintain steady state. 
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Figure 1. 2D slice of the simulation for 600 3 box at 15 Bondi times. Colour represents the entropy, and arrows represent the magnetic 
field vector. The right panel is the equatorial plane (yz), while the left panel a perpendicular slice (xy). White circles represent the Bondi 
radius (rg = 1000). The fluid is slowly moving, in a state of magnetically frustrated convection. A movie of this flow is available in the 
supporting information section of the electronic edition. 



most of our computational effort to the accretion region of 
interest, while also pushing the (periodic) outer boundary 
conditions far away from this region. The inner expansion 
factor Si n is therefore an important numerical parameter, re- 
lated to both the grid's resolution and its anisotropy where 
we care most about the flow. 

Within our inner accretion region, magnetic fields are 
reconnected (relaxed to the vacuum solution consistent with 
the external field, see appendix [B| and mass and heat are 
adjusted (invariably, removed) so that the sound speed and 
Alfven velocity both match the Keplerian velocity at Rb- 
The accretion zone is a cube, whose width we hold fixed 
at 15 in units of the local (uniform) grid separation, so we 
define R[ n — 7.5 dx m i n (but note, the volume of this region 
is equivalent to a sphere of radius 9.3dx m i n .) We consider it 
too costly to vary the numerical parameter i?i n /da; m i n . 

Our grid geometry imposes a local dimensionless reso- 
lution parameter 



max; (dx% 



(4) 



(the maximum being over coordinate directions), which de- 
pends both on radius and on angle within the simulation 
volume. At the inner boundary 5R ~ 7.5 — 9.3; 5ft increases 



to nearly s in 



10 at Rb, then decreases toward s t in 



the exterior region. In Sj3]we report the effective resolution 
at the Bondi radius, 5Rs = ^R(Rb), along with our results. 



2.3 Computational implementation 

Our simulations were performed on the Canadian Institute 
for Theoretical Astrophysics Sunnyvale cluster: 200 Dell 
PE1950 compute nodes; each node contains 2 quad core In- 



tel(R) Xeon(R) E5310 @ 1.60GHz proces sors, 4GB of RA M, 
and 2 gigE network interfaces. The code (|Pen et al.ll2003T ) is 
a second-order accurate (in space and time) high-resolution 
total variation diminishing (TVD) MHD parallelized code. 
Kinetic, thermal, and magnetic energy are conserved and 
divergent of magnetic field was kept to zero by flux con- 
strained transport. There is no explicit magnetic and vis- 
cous dissipation in th e code except on the grid scale. We 
used an MPI version l|Kaeppeli et al.1 120091 ). and 216 CPU 
cores were used to compute a 300 3 box, using OpenMP with 
8 cores per node. The 600 3 simulation was performed on the 
SciNet cluster using 1000 cores over 125 nodes 0. 

A parallel effort to implement the code on economic 
graphics processing units (GPUs) is in progress , which will 
allow larger and longer simulations in the future (|Pang et al.l 
l2010h . 



3 SIMULATIONS AND RESULTS 

Our suite of simulations is described in Table [T] along with 
some selected results. We independently varied the magne- 
tization, rotation, and dynamic range of the flow, as well 
as the effective resolution at Rb- In order to suppress the 
lingering effects of our initial conditions, we ran each simu- 
lation for long enough that a total mass equivalent to all the 
matter initially within Rb was eventually accreted, before 
assessing the flow structure. Because most of our runs ex- 
hibited a significant suppression of the mass accretion rate 
M relative to the Bondi value, this constraint required us 



http://www.scinet.utoronto.ca/ 
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Figure 2. Density versus radius. The dotted line represents the 
density profile for the Bondi solution, which is the steepest plau- 
sible slope at k = 1.5. The dashed line represents the density 
scaling for CDAF solution, which is the shallowest proposed slope 
with k = 0.5. The solid line is the density profile from one of our 
simulations, which is intermediate to the two. 

to simulate for many ts (typically 20 ts)- This requirement 
put strenuous constraints on our simulations (each of which 
required ~ 3 weeks to complete), and will be a serious limita- 
tion on any future simulations performed at higher dynamic 
range. 




Figure 3. log(/3), entropy and radial velocity versus radius. The 
dashed line v r /c s represents the radial velocity in units of mach 
number. The dots v r /c ms represent the radial velocity in units of 
magnetosonic mach number. The solid line is the entropy, and we 
see the entropy inversion which leads to the slow, magnetically 
frustrated convection. Inside the inner boundary, the sound speed 
is lowered, leading to the lower entropy. The + symbols are the 
magnetic field strength, 0. 



version to Bondi inflow if magnetic fields are suddenly elim- 
inated; we have not repeated this experiment.) 



3.1 Character of saturated accretion flows 



3.1.1 Lack of rotational support 



Figure [T] shows the 2D slices for the simulation of our high- 
est resolution 600 3 box at 15 Bondi times (case 25) 0- The 
remaining Figures are all based on case 10, which is most 
representative of the whole set of simulations. Figures[5]and 
[3] display the spherically-averaged properties, figure [5] shows 
the spherically-averaged density of the run; figure [3] shows 
the spherically-averaged radial velocity, /3 and entropy (nor- 
malized to the Bondi entropy). The entropy inversion is 
clearly visible, which leads to the slow, magnetically frus- 
trated convection. 

We draw several general conclusions from the runs listed 
in Table [T] 

- In the presence of magnetic fields, the flow develops a 
super-adiabatic temperature gradient and flattens to k ~ 1. 
Gas pressure remains the dominant source of support at all 
radii, although magnetic forces are always significant at the 
inner radius. 

- Mass accretion diminishes with increasing dynamic 
range, taking values M ~ (2 - 4)Af s (i? in /i? s ) 3/2 " fc . 

- Even significant rotation at the Bondi radius has only 
a minor impact on the mass accretion rate, as the flows do 
not develop rotationally supported inner regions. 

- Our results depend only weakly on the effective resolu- 
tion Kb. 

- In the absence of magnetic fields and rotation, a Bondi 



flow develops. (|Pen et al.l 20031 further demonstrated a re- 



3 Movies are also available in various formats at 
http://www.cita.utoronto.ca/~pen/MFAF/blackhole_movie/index 



The non-rotating character of the flow casts some doubt on 
models which depend on equatorial inflow and axial out- 
flow. Our nonrelativistic simulations cannot rule out an ax- 
ial outflow from a spinning black hole, but they certainly 
show no tendency to develop rotational support in their in- 
ner regions, even after many tens of dynamical times. In a 
rotating run, angular momentum is important at first, in 
preventing the accretion of matter from the equator. Axial, 
low-j material does accrete, but some of it s hocks and drives 
an outflow along the equator (as reported bv lPen et al 1 l2003l 
and lProea fc Be gclman 2003). After a few ts this quadrupo- 
lar flow disappears, leaving behind the nearly hydrostatic, 
slowly rotating envelope which will persists for our entire 
simulation time, i.e. tens oUb- We attribute the persistence 
of this rotational profile to magnetic braking, as the Alfven 
crossing time of the envelope is always shorter than its ac- 
cretion time. Magnetic fields thus play a role here which 
is rather different than in simulations which start from a 
rotating torus, where the magneto-rotational instability is 
the controlling phenomenon; the critical distinction is the 
presence of low-angular-momentum gas. 

Unlike compact object disks, which accrete high- 
angular-momentum material and are guaranteed to cool in 
a fraction of their viscous time, the GCBH feeds upon low 
low-angular-momentum matter, and its accretion envelope 
cannot cool. For both of these reasons it is not surprising to 
discover a thick, slowly rotating accretion envelope rather 
than a thin accretion disk. We stress that global simulations, 
which resolve the Bondi radius and beyond and continue for 
many dynamical times, are required to capture the physical 
htmprocesses which determine the nature of the flow. 
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Run ^ l+s in 5R S ft |f ^ 



1 


500 


67 


1.023 


40.15 


00 





8 


1.02 


1.5047 


2 


250 


33 


1.013 


59.29 


00 





3 


1.10 


1.5273 


3 


125 


17 


1.013 


48.11 


100 





6-20 


0.49 


1.2482 


4 


250 


33 


1.013 


59.29 


100 





6-20 


0.31 


1.1650 


5 


500 


67 


1.023 


40.15 


100 





6-20 


0.22 


1.1399 


6 


1000 


133 


1.0315 


30.82 


100 





6-10 


0.16 


1.1253 


7 


250 


33 


1.013 


59.29 


1 





6-20 


0.15 


0.9574 


8 


250 


33 


1.013 


59.29 


10 





6-20 


0.26 


1.1147 


9 


250 


33 


1.013 


59.29 


1000 





6-20 


0.40 


1.2379 


10 


250 


33 


1.013 


59.29 


100 


0.1 


6-20 


0.289 


1.1450 


11 


250 


33 


1.013 


59.29 


100 


0.5 


6-20 


0.286 


1.1420 


12 


250 


33 


1.013 


59.29 


100 


1.0 


6-20 


0.31 


1.1650 


13^ 


62.5 


33 


1.06 


14.24 


100 





6-20 


0.30 


1.1557 


14c 


125 


33 


1.037 


28.94 


100 





6-20 


0.33 


1.1829 


15~ 


250 


33 


1.013 


59.29 


00 


0.1 


6-20 


0.615 


1.3610 


16 


250 


33 


1.013 


59.29 


00 


0.5 


6-20 


0.621 


1.3637 


17 


250 


33 


1.013 


59.29 


00 


1.0 


6-20 


0.759 


1.4211 


18 


250 


33 


1.013 


59.29 


1000 


0.1 


6-20 


0.400 


1.2379 


19£ 


250 


33 


1.013 


59.29 


1000 


0.1 


6-20 


0.469 


1.2835 


20" 


250 


33 


1.013 


59.29 


100 


0.1 


6-20 


0.300 


1.1557 


21 


250 


33 


1.013 


59.29 


10 


0.1 


6-20 


0.233 


1.0834 


22 


250 


33 


1.013 


59.29 


1 


0.1 


6-20 


0.188 


1.0220 


23 


250 


33 


1.013 


59.29 


100 





6-20 


0.340 


1.1915 


24^ 
25£ 


500 


67 


1.0315 


31.65 


100 


0.1 


6-20 


0.18 


1.2434 


1000 


58.9 


1.015 


64 


100 


0.1 


6-20 


0.19 


1.0925 



a Values are taken from Equation [5] 
^ case of 75 3 grid resolution 
c case of 150 3 grid resolution 

d 19-23, B field is along [0 1] axis 
e 24-25, B field is along [1 2 0] axis 

J case of 600 3 grid resolution 

Table 1. Simulations described in this paper. Columns: Run number; Maximum resolution relative to the Bondi radius; Radial dynamic 
range within Rb\ grid expansion factor within Rg; effective resolution at Rb\ magnetization parameter; rotation parameter; range of 
simulation times over which flow properties were measured; mean mass accretion rate over this period; and typical density power law 
slope (p oc r~ k ) over this period. 



3.1.2 Dependence on parameters; Richardson 
extrapolation 

We now investigate whether our results for the accretion rate 
can be distilled into a single, approximate expression. It is 
clear from the results in Table Q] that rotation affects the 
accretion rate in a non-monotonic fashion. However as we 
have just noted that rotation plays a minor role in our final 
results, we are justified in fitting only the non-rotating runs. 
Rather than M/Afaondi we fit an effective density slope fc e ff 
defined by 

M _ /R in \ 3 / 2 ~ k ^ 
Msondi \RbJ ' 

There are four major variables: the magnitute of the ambient 
magnetic field (A)), the radial dynamical range (Rb/Rih), 
the resolution of the Bondi scale (5Rs). Our fit is 

fc cff = 1.50 - 0.56/? - 098 + 6.51 ~^ - 0.11$R S 048 ; (6) 

all seven numerical coefficients and exponents were opti- 
mized against the 25 runs in Table [I]. The form of equation 
(0 is significantly better than others we tested, including 



those involving log(i?s/i?i n ) and log(SRs). It predicts the 
entries in Table [T]to within a root-mean-square error of only 
0.017. 

Somewhat unexpectedly, this nonlinear fit to our simu- 
lation output recovers the Bondi solution in the continuum, 
unmagnetized limit (fc c g — > 3/2 as ft — > 00, Rs/Rin — > 00, 
and 5Rb — ► 00 )• Moreover the form of the expression al- 
lows us to extrapolate, in the manner of Richardson ex- 
trapolation, to conditions we expect are relevant to Sgr 
A*: K s ~ 00, R B /R in ~ 10 5 , and ft ~ 1 - 5: then, 
fccff ~ 0.94 - 1.0. 

It is encouraging that this result lies in the vicinity 
of observational constraints, lending additional credence to 
the notion that Sgr A* is surrounded by a "magnetically- 
frustrated" accretion flow. We must recall, however, that 
this is only an extrapolation based on simulations which 
lack potentially important physics such as a relativistic in- 
ner boundary and a non-ideal plasma. The absence of an 
imposed outward convective luminosity is likely to be the 
essential element which allows for a lower value of fc. 
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4 ROTATION MEASURE 

The magnitude of RM constrains the density of the inner 
accretion flow, thereby also constraining the mass accretion 
rate and power law index k. Future observations should pro- 
vide time series of RM(t), a rich data set which encodes im- 
portant additional information about the nature of the flow. 
Our goal in this section will be to characterize RM variabil- 
ity within our own simulations sufficiently well to distinguish 
them from other proposed flow classes. 

We pause first to consider why RM should vary at 
all. The rotat ion of polarization is determined by an inte- 
gral (eq. IA1I rShchcrbakovl 120081 ) which is proportional to 
f n e ~B • dl integrated over the zone of nonrelativistic elec- 
trons. The integral is typically dominated by conditions 
at Rrei, the radius where kT e = m e c 2 . Even if n e is rea- 
sonably constant, B likely will change in magnitude and 
direction as the flow evolves. Given that the dynamical 
time at R re \ is under a day, any strongly convective flow 
should exhibit significa nt day-to-day fluctuations in RM; 
measurements bv lMarrone et ail (|2007h appear to rule this 
out. Rotational support also implies rapid RM fluctuations 
unless B is axisymmetric. In the highly subsonic flow of 
magnetically-frustrated convection, however, RM may vary 
on much longer time scales. 

Two proposals have been advanced in which RM(i) 
would be roughl y constant. Within th eir simulations of thick 
accretion disks, Sharma et al.l (|2007T l show that trapping of 
poloidal flux lines leads to a rather steady value of RM for 
observers whose lin es of sight are out of the disk plane. 
Sh arma et "all (|2008l ) point to the constancy of RM in the 
steady, radial magnetic configuration which develops due to 
the saturation of the magneto-thermal instability (in the 
presence of anisotropic electron conduction). We suspect 
that noise at the dynamical frequency is to be expected in 
both these scenarios, which need not exist in a magnetically 
frustrated flow. We also note that both scenarios lead to 
systematically low values of RM for a given accretion rate, 
and therefore imply somewhat higher densities than we in- 
ferred from a spherical model; this may be observationally 
testable. 

Our calculation of RM(t) is based on case 10 in Table 
[1] In Figure [4] we plot RM(i) against an analytical estimate 
of its magnitude. For this purpose we estimate RM as, 



RM: 



2irm e 



n e Bdr, 



(7) 



integrated along radial rays (two per coordinate axis) 
through the simulation volume. We neglect the difference 
between this expression and one which accounts for the rel- 
ativistic nature of electrons within R re \. We therefore nor- 
malize RM to the estimate RM es t as, 



RM C 



2c 4 m! 



GMR Ic in e n e (R r , 



3-1V2 



Hfl- 



given by equation (|A5|) with F(k,kr) — » 1, (cos(0)) — > 1/2, 
P — > 10, and k — ¥ 1. Because we do not calculate elec- 
tron temperature within our simulations, we have the free- 
dom to vary R rc \ and to probe the dependence of coher- 
ence time on this parameter. In practice we chose R IC \ = 
(17, 26, 34, 43)<fa m in in order to separate this radius from 
the accretion zone (7.55x m i n ) and Bondi radius (250(5x m i n , 




Figure 4. Rotation measure vs time (in units of tg). We chose 
^rd = 17, corresponding to R TC i/R^ =0.068. Six lines represent 
three axes: upper set is X (centered at +3), center is Y (centered 
at 0) and lower is Z (centered at -3), with positive and negative 
directions drawn as solid and dashed lines, respectively. 



J 




Figure 5. PDF of RM in Figure[4] The dashed line represents a 
Gaussian distribution. The horizontal axis has been normalized 
by the standard deviation in figure [4] ctrm = 0-63. 



in this case). Figure |4] illustrates RM(t) along each coordi- 
nate axis with the case for R IC \ — 17<5:Tmin). As this figure 
shows, RM changes slowly and its amplitude agrees with our 
estimate RM cs t. In our simulations, we can measure the full 
PDF, shown in figure [5] 

We can ask how well a single measurement of RM con- 
strains the characteristic RM, say the ensemble-averaged 
root-mean-square value RM rms . This is a question of how 
well a standard deviation is measured from a single obser- 
vation. From figure [5] we see that the distribution in our 
simulations is roughly Gaussian with standard deviation 
(7rm = 0.63RM cst . One needs to apply Bayes' Theorem to 
infer the variance of a Gaussian from N independent mea- 
surements: 



/ 2 \ 1/2 
ARM rms = ( — J o rm 



(9) 



To date, no sign change in RM has been observed, suggest- 
ing that we only have one independent measurement. Esti- 
mating RM rms from a single data point requires a Bayesian 
inversion. Estimating from our simulation with a flat prior, 



8 Bijia PANG et al. 




Figure 6. The rotation measure integrant pB r vs radius and time. The central dark bar represents the inner boundary, the vertical axis 
is the Z axis. The horizontal axis is time, in units of ts; Greyscale represents sign(B r ) y p\B r \, which was scaled to be more visually 
accessible. The coherence time is longer at large radii and at late times. Several Bondi times are needed to achieve the steady state 
regime. 



the 95% confidence interval for the ensemble characteristic 
RM given the one data point spans two orders of magnitude! 

In other words, if in fact RM rms = 5.4 x 10 6 , it is not 
very surprising that we have observed RM~ —5.4 x 10 5 . 
The maximum likelihood estimate is RM rms = RM. The 
95% upper bound is RM rms = 33RM, and the lower bound 
is RM rms = 0.33 RM. More data is essential to constrain 
this very large uncertainty. 

A visual description of the RM integrand through the 
flow is shown in figure [6] The time variability time scale is 
shorter at small radii, and shorter at the beginning of the 
simulation. Simulations of many Bondi times with bound- 
aries many Bondi radii away are necessary to see the char- 
acteristic flow patterns. 

To be more quantitative, we plot in Figure [7] the auto- 
correlation of RM(t) for different i? rc i. We define the coher- 
ence time r to be the lag at which the autocorrelation of 
RM falls to 0.5. 

The actual RM radius R Ic i is not resolved in our sim- 
ulations. In order to extrapolate to physically interesting 
regimes, we fit a trend to our limited dynamic range. The 
characteristic variability time scale is given by the flow 
speed, so t a R^ cl p(R Te i) / M . For our characteristic values 
k ~ 1, we have r oc R^ cl , which we fit to our coherence time, 
shown in figure [8] 

For density profiles shallower than Bondi, the charac- 
teristic RM time scale r is significantly longer than the dy- 
namical time (r ~ (f?rei/^B) 3 ^ 2 ts). In our fit, it is given by 
the accretion time 









dot is Rin=43 




-O v dash is Rin=34 




N * s v "~ = * s dash dot is Rin=26 




— ^\ \0 \ solid is Rin=1 7 











i°g,„(t, /t D ) 

M 1tr lags B' 

Figure 7. Autocorrelation for Figure [4] X axis represents 
time lags; Y axis represents autocorrelation for different -Ri n . 
The dotted, dashed, dashed-dot and solid lines correspond to 
R in = 43, 34, 26, 17 respectively. 



(10) 



with a relatively large dimensionless prefactor. This indi- 
cates a coherence time of order one year for the conditions 
at Sgr A*. The actual value of R TC \ is uncertain by a factor 
of 6, so the expected range could be two months to a year. 

This is sufficiently distinct from one day that the dis- 
tinction between frustrated and dynamical flows should be 
readily apparent, once observations span year long baselines. 
We will discuss this point more in Section f5. II below. 
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Figure 8. RM coherence time r as a function of the inner trun- 
cation radius i? re i; points refer to R TB \ = 17, 26, 34 and 43. The 
bootstrap error of 0.17 dex is based on the six data, two for 
each coordinate direction, at each ii rc i. The normalization for 
-Rrd = Rb is log 10 (ti ags /i s ) = 2.15. 

5 DISCUSSION 

In this section we wish to revisit several of the physical pro- 
cesses which are missing from the current numerical simula- 
tions: stellar winds from within Rb, the transport of energy 
and momentum by nearly collisionless electrons, and the in- 
ner boundary conditions imposed by a central black hole. 

Stellar wind input. Our simulations account for the ac- 
cretion of matter from outside the Bondi radius inferred 
from X-ray observations, but not for the direct input of mat- 
ter from in dividual stars in the vicinity of the black hole. 
lLoebl I 2004) raises the possibility that individual stars may 
in fact dominate the accretion flow. The wind from a sin- 
gle star at radius r dominates the flow when its momentum 
output M w v w satisfies 

M w v w > 4Tvr 2 p{r) 3.3(1O~ 5 M yr" 1 )(1000kms _1 ) (11) 

where the evaluation is for a model consistent with RM con- 
straints, in which density follows nu — 10 7 ' 3 (r/i?s) _1 cm _ ' ! 
and pressure follows p ~ 10 4 (r/7is) -2 dyncm -2 - note that 
the criterion is independent of radius for k = 1. The re- 
quired momentum output, equivalent to 10 6 ' 2 L©/c, is well 
above the wind force of any of the OB stars observed within 
Rb- While stars within Rb add fresh matter faster than it 
is accreted by the hole, we can be confident that no single 
star dominates the flow. If the density slope is substantially 
more shallow, for example in a CDAF with k — 1/2, the 
steller winds would be a more important factor. 

Collisionless transport. In the context of a dilute plasma 
where Coulomb collisions are rare, electron thermal con- 
duction has the potential to profoundly alter the flow pro- 
file. The importance of this effect depends on the electrons' 
ability to freely str eam down their temperature gradient 
l|Sharma et al.l [200cf) . despite the wandering and mirroring 
induced by an inhomogeneous magnetic field. The field must 
be weak for the magneto-thermal instability to develop, yet 
weak fields are less resistive to tangling. The thermal con- 
duction is expected to be strongest in the deep interior of 
the flow. If electrons actually free stream inside of 1000 
Schwarzschild radii, the electrons could be non-relativistic 
all the way to the emission region, changing the interpreta- 



tion of the RM. This would favour even shallower density 
profiles, for example the CDAF models. In such a model, 
the RM might be expected to vary on time scales of min- 
utes, which appears inconsistent with current data. If, on 
the other hand, the free streaming length is short on the in- 
side, it more likely places the fluid in an ideal regime for the 
range of radii in our simulations. We therefore remain ag- 
nostic as to the role of thermal conduction in hot accretion 
flows, although it remains a primary caveat of the current 
study. Observations of time variability of RM will substan- 
tially improve our unstanding. 

Black hole inner boundary. Our current inner boundary 
conditions do not resemble a black hole very closely, apart 
from the fact that they also allow gas to accrete. As the inner 
region dominates the energetics of the flow, we consider it 
critical to learn how the black hole modifies our results. We 
are currently engaged in a follow-up study with a relativistic 
inner boundary, to be described in a future paper. 

5.1 Observational Probes 

RM can be measured by several techniques. Currently, ef- 
forts have concentrated at high frequencies, v ~ 200 — 300 
GHz (jMarrone et alj l2006t ). where the polarization angle 
varies slowly with frequency. Accurate measurements over 
long time baselines allows discrimination between models. 
At high frequencies, the SMA and ALMA would allow a 
steady synaptic monitoring program. The full time 2 point 
correlation function extends the measurement space by one 
dimension. 

At lower frequencies, high spectral resolution is needed 
to resolve the winding rate, which is now tractable with 
broad band high resolution instruments such as the EVLA 
and ATCA/CABB. The higher winding rate would allow a 
much more sensitive measurement of small changes in the 
RM, which would also be a descriminant between models. 
The challenge here is that the polarization fraction drops 
signficantly with frequency, requiring a more accurate in- 
strumental polarization model. On the other hand, the very 
characteristic A 2 dependence of RM should allow a robust 
rejection against instrumental effects. 

At lower frequency, the spatial extent of the emission re- 
gion is also expected to increase. When the emission region 
approaches the rotation measure screen, one expects depo- 
larization. Direct polarized VLBI imaging could shed light 
on this matter. This is complicated by interstellar scattering, 
which also increases the apparent angular size. The changing 
emission location a s a function of frequency may complicate 
the R M inferences dFish et al . 2009; lBroderick fc McKinnevI 
|2010|) This can skew the actual value of the infered RM, re- 
sulting in an underestimate. The sign of RM would generi- 
cally be a more robust quantity, and looking for changes in 
the sign of RM could be a proxy for the correlation function. 

A separate approach is to use other polarized sources 
as probes of the flow. One candidate population is pulsars. 
At the galactic center, interstellar scattering (Lazio 2000) 
smears out the pulses, making them difficult to detect di- 
rectly. But the pulse averaged flux should still be present. 
Over the orbit around the black hole, one can measure the 
time variation of the RM, leading to a probe of the spatial 
RM variations in the accretion flow. Some pulars, such as the 
crab, exhibit giant pulses, which could still be visible despite 
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a scattering delay smearing. These could be used to measure 
the dispersion measure (DM) along the orbit. The GMRT 
at 610 MHz would have optimal sensitivity to detecting the 
non-pulsating emission from pulsars, and be able to decon- 
fuse them from the domin ant synchrotron emission us ing 
rotation measure synthesis ()Brentiens fc de Bruvnll2005l ). 



Fund. The computations for the movies were performed on 
the GPC supercomputer at the SciNet HPC Consortium. 
SciNet is funded by: the Canada Foundation for Innovation 
under the auspices of Compute Canada; the Government of 
Ontario; Ontario Research Fund - Research Excellence; and 
the University of Toronto. 



6 CONCLUSION 

A series of new, large dynamical range secular MHD simula- 
tion are presented for the understanding of the low luminos- 
ity of the supermassive black hole in the Galactic Center. 
These are the first global 3-D MHD simulations which do 
not face boundary conditions at outer radii, and impose in- 
going boundaries at the interior, running for many Bondi 
times. We confirm a class of magnetically frustrated accre- 
tion flows, whose bulk properties are independent of physi- 
cal and numerical parameters, including resolution, rotation, 
and magnetic fields. No significant rotational support nor 
outward flow is observed in our simulations. An extrapola- 
tion formula is proposed and the accretion rate is consistent 
with observational data. 

A promising probe for the nature of the accretion flow 
is the rotation measure, and its time variability. In this com- 
parison, the dominant free parameter is the electron temper- 
ature. We argued that over the plausible range, from thermal 
to adiabatic, this radius varies from 40 to 250 Schwarzschild 
radii. The RM variations in the simulations are intermit- 
tent, requiring many measurements to determine this last 
free parameter. 

We propose that temporal rotation measure variations 
are a generic prediction to distinguish between the wide va- 
riety of theoretical models currently under consideration, 
ranging from CDAF through ADIOS to ADAF. RM is domi- 
nated by the radius at which electrons turn relativistic, when 
the flow is still very subrelativistic, and is thus much fur- 
ther out than the Schwarzschild radius. Most models, other 
than the ones found in our simulations, involve rapidly flow- 
ing plasmas, with Mach numbers near unity. These generi- 
cally result in rapid RM variations on time scales of hours 
to weeks (or in special cases, it can be infinite). In con- 
trast, our simulations predict variability on time scales of 
weeks to years. A major uncertainty in this prediction is the 
poor statistical measure of the standard deviation of RM 
measurement, which requires long term RM monitoring to 
quantify. 

Future observations of RM time variability, or spatially 
resolved measurements using pulsars, will provide valuable 
information. 
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APPENDIX A: ROTATION MEASURE 
CONSTRAINT ON ACCRETION FLOW 

In traversing the accretion flow, linearly polarized radio 
waves of wavelength A are rotated by RM A 2 radians, where 



RM = 



2nm 2 c 4 



n e f(T e )Bcos(6)dl. 



(Al) 



Here f(T e ) is a ratio of modified Bess el functions: 
f(T e ) = K (m e c 2 /k B T e )/K 2 (m e c 2 /k B T e ) JShcherbakovl 
2008), which suppresses RM by a factor oc T~ 2 wherever 
electrons are relativistic. The integral here covers the en- 
tire path from source to observer; is the angle between 
B and the line of sight. This expression is appropriate 
for the frequencies at which RM has been observed; at 
lower frequencies, where pro pagation is "superadiabatic" 
l|Broderick fc Blandfordll2010] ) cos(6>) ->■ ±| cos(0)|. 

We adopt a power-law solution with negligible rota- 
tional support in which p oc r~ h , and the total pressure 
P oc r~ kp with kp = k + 1; moreover we take T e oc r~ kT for 
the relativistic electrons. The hydrostatic equation dP/dr — 
—GMp/r 2 becomes 



P = Pa + P B = 



GM p 
(k + 1) r' 



and with P g = /3P B = (3B 2 /(8tt) 
1.2m p is the mass per electron), 



B 



8tt 



GMp, e n e 



(/? + !)(&+ 1) 



1/2 



(A2) 

p — n e p e (where fi e = 
(A3) 



So long as k > 1/3 (so that RM converges at large radii) 
and k < (1 + 4fcr)/3 (so it converges inward as well), the 
RM integral is set around R vc \. Taking a radial line of sight 
(dl — > dr), we write 




n e f(T e )Bcos(6)dr = F(k,k T ) / n c Bcos(6)dr (A4) 

1 -'-Rrcl 

2 'cos(6))F(k,k T )[n e Br 



3k - 1 

where (cos(#)} encapsulates the difference between the true 
integral what it would have been if = all along the path, 
and F{k,kr) encapsulates the difference between a smooth 
cutoff and a sharp one. We plot F(k, kr) in Figure I All it 
is of order unity except as kr approaches (3k — l)/4. All 
together, 



RM : 



4e 2 GM (cos(6)} F(k,k 7 
m 2 c 5 3k- 1 



/x e n e (i? re i) 3 



Rr, 



jr(jfe + 1)08 + 1) R s 



1/2 



.(A5) 



Figure Al. The logarithm of the relativistic RM factor, 
log 10 F(k, kx)- The true RM integral is modified by a factor 
F(k,kx) relative to an estimate in which the nonrelativistic for- 
mula is used, but the inner bound of integration is set to the 
radius R Ic i at which electrons become relativistic; see equation 
[All 



To estimate n e (i? rc i) from RM, one must make assump- 
tions about the uncertain parameters /3, (cos(#)), kr, and 
R TE i/Rs; then k can be derived self-consistently from ob- 
servations n e (R B ) and RM. Our fiducial values of these pa- 
rameters are 10, 0.5, 0.5 and 100, respectively, of which we 
consider the last to be the most uncertain. We now discuss 
each in turn. 

Although the magnetization parameter j3 could conceiv- 
ably take a very wide range of values, we consistently find 
/3 ~ 10 in our simulations, with some tendency for /3 to 
decrease inward. We consider it unlikely for the flow to be 
much less magnetized, given the magnetization of the galac- 
tic center and the fact that weak fields are enhanced in most 
of the flow models under consideration. 

If B wanders little in the region where the integrand is 
large (a zone of width ~ _R rc i around R TC i), and is randomly 
oriented relative to the line of sight (cos(#)) ~ cos(8(R TC \)), 
typically 1/2 in absolute value. If the field were purely 
radial, (cos(#)) would be unity. Conversely if B reverses 
frequently in this region (the number of reversals N r is 
large) then (cos(0)) rma ~ 1/(2 V / A r r + 1) will be small. How- 
ever, N r cannot be too large, or magnetic forces are un- 
balanced. We gauge its maximum value by equating the 
square of the buoyant growth rate, M 2 = [(3-2k)/5]GM/r 3 , 
against the square of the Alfven frequency N 2 v\/r 2 . Not- 
ing that v\ = GM/[2(/3 + l)(k + l)r], we find N 2 ~ 
(2/5)(/3 + l)(fc + 1)(3 - 2k). For /3 = 10 and k = 1 this 
implies (cos(6 , )) rma ~ 0.25: a very minor suppression. We 
can therefore be confident that (cos(#)) = 0.5 to within a 
factor of 2, unless /? 3> 10 for some reason. 

The precise value of kr is not important unless it ap- 
proaches or falls below the minimum value (3fc — 1)/4. If elec- 
tron conduction is very strong this is unavoidable, as rapid 
transport implies kr — 0; however in this case the relativis- 
tic region disappears, as discussed below. Alternately, if rel- 
ativistic electrons are trapped and adiabatic, T e oc p 1 ' 3 and 
kr = k/3; however kr < (3k — l)/4 then requires k < 3/5, 
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which can only be realized within the CDAF model. Finally, 
if electrons remain strongly coupled to ions, kr = 1 and we 
only require k < 5/3. 

The location at which electrons become relativistic, 
- Rrei, i s quite uncertain. Models such as those of I Yuan et alj 
(2003), in which electrons are heated while advecting inward, 
predict R Tcl ~ 10 2 R s . The maximum conceivable _R re i cor- 
responds to adiabatic compression of the electrons, inward 
from the radius at which they decouple from ions; this yields 
about 500/(1 + k)Rs . If conduction is very strong, however, 
electrons should remain cold throughout the flow; in this 
case we should replace R Ic i/Rs — > 1 and F(k,kr) — > 1 in 
equation (|A5|) . 

Adopting our fiducial values for the other variables, and 
taking F(k, kr) — > 1 for lack of knowledge regarding kr, we 
may solve for the self-consistent value of k which connects 
the density at Rb with n e (-R r el) derived from equation IA5I 
We find k -> (0.90, 1.23, 1.32) for R Ici /R s -> (200, 100, 1), 
respectively. As noted in the text, the current small set of 
RM measurements allows a two order of magnitude range in 
RM cst , and k ~ 1 is consistent with data. Longer observa- 
tions of time and amplitude will improve the constraints. 



APPENDIX B: 
CONDITIONS 



INNER BOUNDARY 



The inner boundary conditions were determined by first 
solving for the vacuum solution of the magnetic field in- 
side the entire inner boundary cube. Then inside the largest 
possible sphere within this cube, matter and energy were 
removed. 

To simplify the programming, we put the entire inner 
boundary region on one node. This meant that the grid had 
to be divided over an odd number of nodes in each Cartesian 
direction. 



Bl Magnetic field 

In order to determine the vacuum magnetic field solution, 
we use the following two Maxwell equations for zero current: 



V • B = , 

V x B = 



(Bl) 
(B2) 



Equation (|B2|) enables us to write B = V(f>, for some scalar 
function <f>. Combining this with (|B1|) we obtain Laplace's 
equation 



V 2 = O , 



(B3) 



which we solve with Neumann boundary conditions (the nor- 
mal derivative n-\7<f> specified) given byB-non the bound- 
ary of the cube. 

Since the MHD code stores the values of B on the left- 
hand cell faces, we must solve for <j> in cell centers and then 
take derivatives to get the value of B on the cell boundary. 
Let the inner boundary cube be of side length N, consisting 
of cells numbered 1, . . . , N in all three directions. In order 
to simplify the problem we set B ■ n — on five of the six 
faces of the cube, and find the contribution to (j> from one 
face at a time. 



Suppose B n = on all of the faces except the i — N+l 
face (i.e., B^ +1 ' 3 ' can be non-zero). The Laplace equation 
(|B3|) with Neumann boundary conditions only has a solution 
if the net flux of field into the cube is zero. Since all of the 
boundary conditions are zero except for the i — N + 1 face, 
that face must have a net flux through it of zero. Defining 



(B4) 



R JV+l _ 1 V^V^ r>N+l,j,k 

a * ~ jv 2 2—< 

3=1 k=l 

to be the average of B x on the i = N + 1 face, and letting 

b N+l,J,k = B N+l,j,k _ ^iV+T b N+l,j,k can be uged ag the 

boundary condition and B^ +1 will be added in later. 

We use separation of variables to solve for <f>. Set <j>^ = 
X i Y j Z h , substitute into (|B3j, and rearrange to get 

X i+1 — 2X i + X^ 1 Y J+1 - 2Y J ' + Y 1 ' 1 
Xi + ~ h 

Now let 

yj+i _ 2Y 3 + Y^ 1 
Y~J 

Z k+i _ 2Z k + Z k-i 

Solving equations (|B6[) and ()B7|I with the boundary condi- 
tions, 



= — rj , and 



(B6) 



(B7) 



Yl = cos 



7B7T(j - i) 



N 



k rwr(fc-i) 
Z„ — cos — 



N 



2 ■ 2 17177 2 2 n7r 

ri m = 4 sin — — , lj„ = 4 sm — — ■. 



(B8) 
(B9) 

Substituting (TB6)) . (TBT)) . and (TB9)) into (TB5)) . and solving, 
yields 

0!mnVr(j — i) 



2N 



where 



2A 

7T 



_ . 717T , m7T 

OLmn = arcsmhw sm — — + sm 



2N 2N 
Finally, putting this all together, 

JV-l N-l 



(B10) 
(Bll) 



m - EE- 4 - cos ™ u ~~ k) cos nAk - |} x 



N 



N 



m— n — 



x cosh 



«mn7r(i — |) 



N 



(B12) 



and define Aoo = 0. 

To determine the coefficients A mn we add in the final 
boundary condition (i = N + 1), and get 

4 1 



A mn — 



N 2 2sinh(a m „7r) sinh(a mn 7r/2iV) 

AT jv /• u 

h N+l,j,k _ m7r U ~ 2> s 



x e E b * ' ' J " c " s ' 

j=l k=l 

rm(k — i) 

x cos — — 

N 



N 



(B13) 



A similar calculation may be performed for the case when 
the i = 1 boundary has non-zero field. After finding the 
contribution from each face, store their sum in 6. 
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(B19) 
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Figure Bl. Vacuum solution of the magnetic field is calculated 
in the central region. The field lines outside of the central region 
show the boundary condition. 



We then set p to O.lp. p and p were assigned minimum values 
of 0.1 times the average value of p outside of the sphere, and 
0.001, respectively, to ensure stability. 

This paper has been typeset from a TpX/ WFfjfi. file prepared 
by the author. 



APPENDIX C: SUPPORTING INFORMATION 

Movie 1. Animation of magnetically frustrated convection 
simulation. 

The qualitative behavious of the accretion flow is best 
illustrated in the form of a movie. This movie shows case 25. 
The raw simulation used 600 3 grid cells. The Bondi radius is 
at 1000 grid units, where one grid unit is the smallest central 
grid spacing. The full box size is 8000 3 grid units. Colour 
represents the entropy, and arrows represent the magnetic 
field vector. The right side shows the equatorial plane (yz). 
the left side shows a perpendicular plane (xy). The moving 
white circles represent the flow of an unmagnetized Bondi 
solution, starting at the Bondi radius. On average, the fluid 
is slowly moving inward, in a state of magnetically frus- 
trated convection. Various other formats can also be seen at 



To deal with the subtracted cube face field averages, let |http://www.cita.utoronto.ca/~p en/MFAF/bl ackhole-movie/index.html 



^ k = Bl> k i + Bl lk j + Bl jl k + 



+ - 



B, 



N+l,j,k 



B. 



ljk 



2N 



+ - 



2N 

and add this to <j>. <j>o is the potential of a cube where each 
face has the uniform magnetic field given by the average of 
the magnetic field on the corresponding face of the inner 
boundary cube. 
To find B, set 

(B15) 
(B16) 
(B17) 



^ X 






fjijk 

JDy — 













In Figure IB1I we used the magnetic field solver with a 
boundary condition consisting of field going in one side and 
out an adjacent side of the box. This boundary condition 
tests both the <f>o component of the solution (since faces have 
non-zero net flux) as well as the Fourier series component 
(since faces have non-constant magnetic field). 



B2 Density and pressure 

Inside the largest possible sphere that can be inscribed 
within the inner boundary cube, we adjust the density and 
pressure so that the Alfven speed and the sound speed are 
both equal to the circular speed. We accomplish this by set- 
ting 



